home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / linalg / householdercomplex.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  4.2 KB  |  142 lines
  1. /* linalg/householdercomplex.c
  2.  * 
  3.  * Copyright (C) 2001 Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include <config.h>
  21. #include <stdlib.h>
  22. #include <gsl/gsl_math.h>
  23. #include <gsl/gsl_vector.h>
  24. #include <gsl/gsl_matrix.h>
  25. #include <gsl/gsl_blas.h>
  26. #include <gsl/gsl_complex_math.h>
  27.  
  28. #include <gsl/gsl_linalg.h>
  29.  
  30. gsl_complex
  31. gsl_linalg_complex_householder_transform (gsl_vector_complex * v)
  32. {
  33.   /* replace v[0:n-1] with a householder vector (v[0:n-1]) and
  34.      coefficient tau that annihilate v[1:n-1] */
  35.  
  36.   const size_t n = v->size ;
  37.   
  38.   if (n == 1)
  39.     {
  40.       gsl_complex alpha = gsl_vector_complex_get (v, 0) ;      
  41.       double absa = gsl_complex_abs (alpha);
  42.       double beta_r = - (GSL_REAL(alpha) >= 0 ? +1 : -1) * absa ;
  43.  
  44.       gsl_complex tau;
  45.       GSL_REAL(tau) = (beta_r - GSL_REAL(alpha)) / beta_r ;
  46.       GSL_IMAG(tau) = - GSL_IMAG(alpha) / beta_r ;
  47.  
  48.       {
  49.         gsl_complex beta = gsl_complex_rect (beta_r, 0.0);
  50.         gsl_vector_complex_set (v, 0, beta) ;
  51.       }
  52.       
  53.       return tau;
  54.     }
  55.   else
  56.     { 
  57.       gsl_complex tau ;
  58.       double beta_r;
  59.  
  60.       gsl_vector_complex_view x = gsl_vector_complex_subvector (v, 1, n - 1) ; 
  61.       gsl_complex alpha = gsl_vector_complex_get (v, 0) ;            
  62.       double absa = gsl_complex_abs (alpha);
  63.       double xnorm = gsl_blas_dznrm2 (&x.vector);
  64.       
  65.       if (xnorm == 0 && GSL_IMAG(alpha) == 0) 
  66.         {
  67.           gsl_complex zero = gsl_complex_rect(0.0, 0.0);
  68.           return zero; /* tau = 0 */
  69.         }
  70.       
  71.       beta_r = - (GSL_REAL(alpha) >= 0 ? +1 : -1) * hypot(absa, xnorm) ;
  72.  
  73.       GSL_REAL(tau) = (beta_r - GSL_REAL(alpha)) / beta_r ;
  74.       GSL_IMAG(tau) = - GSL_IMAG(alpha) / beta_r ;
  75.  
  76.       {
  77.         gsl_complex amb = gsl_complex_sub_real(alpha, beta_r);
  78.         gsl_complex s = gsl_complex_inverse(amb);
  79.         gsl_blas_zscal (s, &x.vector);
  80.       }
  81.       
  82.       {
  83.         gsl_complex beta = gsl_complex_rect (beta_r, 0.0);
  84.         gsl_vector_complex_set (v, 0, beta) ;
  85.       }
  86.       
  87.       return tau;
  88.     }
  89. }
  90.  
  91. int
  92. gsl_linalg_complex_householder_hm (gsl_complex tau, const gsl_vector_complex * v, gsl_matrix_complex * A)
  93. {
  94.   /* applies a householder transformation v,tau to matrix m */
  95.  
  96.   size_t i, j;
  97.  
  98.   if (GSL_REAL(tau) == 0.0 && GSL_IMAG(tau) == 0.0)
  99.     {
  100.       return GSL_SUCCESS;
  101.     }
  102.  
  103.   /* w = (v' A)^T */
  104.  
  105.   for (j = 0; j < A->size2; j++)
  106.     {
  107.       gsl_complex tauwj;
  108.       gsl_complex wj = gsl_matrix_complex_get(A,0,j);  
  109.  
  110.       for (i = 1; i < A->size1; i++)  /* note, computed for v(0) = 1 above */
  111.         {
  112.           gsl_complex Aij = gsl_matrix_complex_get(A,i,j);
  113.           gsl_complex vi = gsl_vector_complex_get(v,i);
  114.           gsl_complex Av = gsl_complex_mul (Aij, gsl_complex_conjugate(vi));
  115.           wj = gsl_complex_add (wj, Av);
  116.         }
  117.  
  118.       tauwj = gsl_complex_mul (tau, wj);
  119.  
  120.       /* A = A - v w^T */
  121.       
  122.       {
  123.         gsl_complex A0j = gsl_matrix_complex_get (A, 0, j);
  124.         gsl_complex Atw = gsl_complex_sub (A0j, tauwj);
  125.         /* store A0j - tau  * wj */
  126.         gsl_matrix_complex_set (A, 0, j, Atw);
  127.       }
  128.       
  129.       for (i = 1; i < A->size1; i++)
  130.         {
  131.           gsl_complex vi = gsl_vector_complex_get (v, i);
  132.           gsl_complex tauvw = gsl_complex_mul(vi, tauwj);
  133.           gsl_complex Aij = gsl_matrix_complex_get (A, i, j);
  134.           gsl_complex Atwv = gsl_complex_sub (Aij, tauvw);
  135.           /* store Aij - tau * vi * wj */
  136.           gsl_matrix_complex_set (A, i, j, Atwv);
  137.         }
  138.     }
  139.       
  140.   return GSL_SUCCESS;
  141. }
  142.